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1. Introduction 

Advances in computer hardware and numerical methods in conjunction with care- 
fully designed computer programs have made meaningful numerical simulation of 
wall-bounded turbulent flows possible. The physical realism of the resultant 
computer-generated data has been validated by detailed structural and statistical 
comparisons with experimental measurements. These calculations have proven to be 
a very useful complement to the laboratory experiments. 

This paper reviews some recent developments in three-dimensional, time- 
dependent numerical simulation of turbulent flows bounded by a wall. We shall be 
considering both direct and large-eddy simulation techniques within the same com- 
putational framework. In the following section, we have outlined the governing 
equations. In Section 3, the computational spatial-grid requirements as dictated 
by the known structure of turbulent boundary layers are presented. Next, we review 
the numerical methods currently in use. Some of the features of these algorithms, 
including spatial differencing, time advancement, and data management, will be dis- 
cussed in some detail. In Section 5 we provide a selection of the results of 
recent calculations of turbulent channel flow, including the effects of system 
rotation and transpiration on the flow structure. Finally, in Section 6 we shall 
make our concluding remarks. 

2. Dynamical Equations 

To date, attention has been largely focused on the incompressible flows 
governed by the Navier-Stokes equations: 
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where Re is the Reynolds number. To satisfy the incompressibility constraint 
(lb), certain numerical techniques use the Poisson equation for the dynamic pressure: 
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obtained from application of the divergence operator to Eq. (la). 

In the direct simulation (DS) approach, the above equations are solved numeri- 
cally with appropriate boundary conditions. Aside from errors due to numerical 
implementation, no further approximations are required. In the large-eddy simula- 
tion (LES) technique, the dependent variables are the resolvable portion of the 
velocity and pressure field. Every flow variable f is decomposed to large-scale 
and subrid-scale components. 


f = f + f’ (2) 

The large-scale component is defined by 

f(x) = y'c(x,x’) f(x’) dx’ (3) 

where G is a filter function with a characteristic length A, which is a function 
of the computational grid resolution. Applying the filtering operation (3) to 
Eqs. (la), (lb), or (Ic) leads to the exact equations for the large-scale field. 

The major difference between the filtered and unfiltered equations of the direct 
simulation is the inclusion of the additional terms associated with the subgrid-scale 
stresses (SGS) in the governing equations for the large eddies. These terms are 
modeled to close the system of equations. Different eddy-viscosity models D, 2, 3D 
as well as multi-equation models CaD have been successfully used to relate SGS 
stresses to the resolvable turbulence. These models should display an important fea- 
ture: as the grid resolution is refined, the characteristic length of the SGS eddies 

becomes smaller. 

In this paper, most of the discussion will be in reference to a Cartesian coor- 
dinate system. The x and z (xi , X3) axes are parallel to the wall, with x 
increasing in the mean-flow direction. The y axis is normal to the wall. We shall 
primarily discuss numerical simulation of flows that are homogeneous or nearly homo- 
geneous in the x and 2 directions. This is the area where most of the current 
effort has been concentrated. Some computations of unidirectional flows in cylindri- 
cal geometries have been performed [2, SD- Currently, however, numerical simulation 
of the turbulent flows in complex geometries involving generalized coordinate systems 
has not been undertaken. Other notions used in this paper include: 5, the channel 
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half-width or boundary-layer thickness; U^, free-stream velocity or centerline 
velocity in channel; shear velocity; U^, average mean velocity; Re, the 
Reynolds number based on and 6; Re.^, the Reynolds number based on u^ 

and 6; v, the kinematic viscosity; y^, the distance to the wall; and N^, 
number of grid points in the x j^-dir action . 

3. Spatial Resolution Requirements 

Generally, numerical simulation of wall-bounded turbulent shear flows requires 
a large number of grid points in all spatial directions. This requirements is much 
more stringent than the corresponding one for free-turbulent shear flows (e.g. , jets 
and wakes). The difference stems from physical observations that locally large 
eddies near the wall are much smaller than those away from the wall. Moreover, in 
free-turbulent flows, large-scale structures exhibit an appreciable degree of 
Reynolds number independence, whereas, near the walls, the important large-scale 
structures decrease in size with increasing Reynolds number. 

In the direction normal to the solid boundary, one should and can easily dis- 
tribute grid points with variable spacings. A sufficient number of grid points 
should be placed near the wall to resolve the viscous sublayer and the buffer layer. 
As the Reynolds number increases, more points are required in this region. The grid 
size can be increased in the regions away from the walls; however, it should be 
bounded by the Prandtl mixing length ('*^0.1 6). 

In the lateral (spanwise) direction, the required number of grid points can be 
prohibitively large. The difficulty arises due to the fine spacing of the streaky 
structures C6]] in the vicinity of the wall. Kline and his co-workers have shown 
that streaks play a significant role in the production and dynamics of turbulence in 
the entire flow field. Therefore, it is important that the numerical simulation of 
wall-bounded turbulent shear flows resolve these structures or account for their 
effects. For a limited range of Reynolds numbers, laboratory observations, as well 
as some quantitative measurements indicate that the mean spacing of the streaks, 
is about 100 wall units, i.e., = A^u^/v 100, and their most probable spacing is 

about 80 wall units. The mean width of the smaller of the high- and low-speed 
streaks can be at most 50 v/u^. In fact, from the measurements of Blackwelder and 
Kaplan , one can deduce that the mean width of the high-speed wall-layer struc- 
tures is about 20-40 wall units. These values are based on an ensemble of measure- 
ments, and, at a given instant, structures with smaller (as well as larger) widths 
are formed. Therefore, in order to capture the wall structures at their proper 
scale, it is not unreasonable CsD to require that the computational grid resolution 
in the spanwise direction be fine enough to resolve eddies with a spanwise extent of 
20 wall units. In the numerical integration of the governing nonlinear equations, if 
we assume that at least four grid points are required to represent an eddy and its 
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evolution for a short period of time, the computational grid size in the z-direction 
should be about 5 wall units, i.e., hf - 5. It should be pointed out that this 
estimate is based on experimental data for moderately low Reynolds number turbulent 
flows and may not apply at very high Reynolds numbers. 

In the absence of physical boundaries in the spanwise direction, the extent of 
the computational domain in this direction, L^, should be large enough that arti- 
ficialities introduced by the use of periodic boundary conditions do not seriously 
affect the statistics of the flow. Based on two-point correlation measurements [91, 

L should be at least three times the boundary-layer thickness [3l . VJith these two 
z 

estimates, the required number of grid points in the z-direction, , is 


^ — - — = 0.6 Re 


‘^z 5 /Re 


Using the universal velocity-distribution law ClOJ > may relate Re to Re^ i 


Re = Re^ In Re^ + 5.0 + 


V 7 here k 0.4, E - 0 for channel and pipe flows, and E = 2.8 for the boundary 
layer. Figure 1 shows the required number of grid points in the spanwise direction 
vs. the Reynolds number for channel flows. 



Figure 1. Grid-point requirements in the spanwise 
direction. 
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The computational grid requirement is not as stringent in the streamwise 
direction as in the spanwise direction. Similar considerations of the physics of 
turbulent boundary layers tend to indicate C3] that, in the streamwise direction, 
the required number of grid points is about half that for the spanwise direction.* 

As an example, for the moderate Reynolds number Re = 10^, 64 nonuniformly 
spaced grid points appear to be sufficient in the direction normal to the wall(s). 
For this case, therefore, about 4 x lO^ mesh points are necessary to resolve the 
energetic turbulent structures. The computational effort required for this calcu- 
lation overly taxes the capabilities of the presently available supercomputers. On 
the other hand, low Reynolds number flows, such as the channel flow experiment of 
Eckelmann D-l] with Re = 2800, is definitely within reach of such computers. For 
this flow, less than half a million grid points are required. 

It is emphasized that the spatial resolution requirements just given are based 
on the experimentally determined **large’* eddy sizes. In the above calculations, one 
must use subgrid scale (SGS) models to represent the small dissipative eddies. It 
is difficult to characterize the size of these very small eddies. As a reference, 
however, we may consider the Kolmogoroff length scale n = (v as their 

typical length scale. For channel flows based on the mean dissipation rate per unit 
mass for the entire flow field, n can be expressed in the wall units as 
T\^ = ^ ^ Reynolds number, Re = 2800, that was just con- 

sidered, n“*" is approximately 2 wall units. Based on the wall value of e, the 
limiting value of is slightly less than 1 (n"*" is exactly equal to 1, if only 

the dissipation due to mean motion is considered) . Clearly, if eddies whose extent 
is about in all the spatial directions are to be resolved, the required number 

of grid points is prohibitively large. Thus, it appears that, with the present com- 
puters, direct numerical simulation of wall-bounded turbulent shear flows is not 
feasible, and calculations should incorporate subgrid scale models to represent 
the dissipative eddies. However, this conclusion has been based on using as 

the characteristic size of the dissipative eddies that must be resolved. The 
question arises as to whether eddies of this size make significant contributions to 
the local dissipation rate in turbulent boundary layers. With the available experi- 
mental data, it is difficult to answer this question conclusively. A rough estimate 
of the dissipation spectra obtained from the one-dimensional energy spectra measure- 
ments of Bakewell and Lumley [lH seems to indicate that these eddies do make appre- 
ciable contributions to the local dissipation rate. Perhaps, the easiest way to 
answer this question is to attempt a direct simulation of a very low Reynolds number 


*This estimate is based on the typical streamwise extent of the wall-layer struc- 

trues C12D . In visual studies CbD , the "lifted" sublayer streaks were observed 
to oscillate. Depending on the wavelength of these oscillations, more points in 
the streamwise direction may be necessary. 
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turbulent channel flow and validate it by detail structural and statistical compari- 
son with the available experimental data. 

The demand for a large number of grid points for resolving the wall-layer 
structures can be significantly reduced by the grid-imbedding technique . One 
can place a large number of grid points in the x and z directions only in the 
vicinity of the walls. Since implicit time advancing will generally be used in con- 
junction with derivatives in the normal direction, this approach is not as conveni- 
ent to implement as one where the same number of grid points in the x and z 
directions are used at all the y-locations. However, careful grid imbedding can 
lead to enormous savings on computer time and storage. If at moderate Reynolds num- 
bers the wall-layer structures are to be resolved at their proper scale, grid imbed- 
ding appears to be the only course of action at present. 

Another method for alleviating the need for a large number of grid points in 
the simulation of wall-bounded shear flows is that of Deardor ff [l] and Schumann C2] . 
In this method the flow in the vicinity of the walls is ignored. The calculations 
are carried out to a point in the logarithmic layer where boundary conditions con- 
sistent with the logarithmic velocity distributions are applied. For high Reynolds 
number flows and certain practical problems, this approach is very promising. With 
considerably less effort than is required to extend the calculations to the wall, 
successful comparison of the mean velocity and turbulent stresses with experimental 
data has been obtained. However, the applicability of these empirical boundary 
conditions to other flow situations has not been established. Moreover, the effect 
of perturbations to these boundary conditions on the computed flow field is not yet 
known. If a two-dimensional, time-dependent "wall function" can indeed be con- 
structed with a sufficient degree of generality, then this type of calculation can 
be of considerable practical value for the numerical simulation of high Reynolds 
number, wall-bounded turbulent flows. The calculations that do extend to the wall 
can serve as a viable testing ground to validate the proposed wall conditions. 

A novel and inexpensive method for evaluating these conditions is described by 
Chapman and Kuhn [15] . They calculated the inner region of a turbulent boundary 
layer by specifying space- and time-dependent boundary conditions at the outer edge 
of the viscous sublayer. Considerable care was exercised in assuring that these 
conditions were consistent with the known dynamics of the near-wall turbulence. The 
appropriate "outer" boundary conditions used in this work can be used as wall condi- 
tions in the large-eddy simulation of the outer region of turbulent boundary layers. 

4. Numerical Methods 

In this section we shall discuss the numerical methods used to solve the three- 
dimensional, time-dependent, Navier-S tokes equations for wall-bounded turbulent 
flows. As a result of modeling the subgrid scale stresses, the dynamical equations 
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in the large-eddy simulation approach are somewhat more complicated than Eqs. (1). 
However, when eddy-viscosity models are used, these equations can be recast into a 
form for which virtually identical numerical methods can be used- Equation (la) can 
be written as: 

SU. CiT) 1 

i = u _ JL + i 

9t i 8 x. R 9x.Bx. ^ ^ 

1 J J 

where represents the nonlinear terms (including the subgrid scale terms) and 

the dependent variables (u^,P) are to be interpreted either as the full velocity and 
pressure field in DS , or their resolvable portions in LES- 


4.1 Spatial Derivatives 


Finite-difference and spectral (pseudospectral) methods have been used to 
approximate the partial derivatives 3/3x^. Deardorff [l] and Schumann L2j used 
second-order finite differences in conjunction with staggered grids Ql 6 j in all 
spatial directions. Moin and Kim C3D evaluated partial derivatives in two of the 
spatial directions (xj , X 3 ) pseudospectrally , whereas, derivatives in the direction 
normal to the walls were evaluated by second-order central-difference formulae. 

Orszag and co-workers Cl7, 182, Moin and Kim Cl92 , Kleiser [202 and Taylor and 
Murdock [2l2 computed all the spatial derivatives by pseudospectral methods. When 
using pseudospectral methods, the dependent variables are expressed as a linear com- 
bination of a set of orthogonal functions. Fourier series is the appropriate repre- 
sentation of the flow field in the directions for which periodic boundary conditions 
are used. In other dierctions, orthogonal polynomial decompositions generally lead 
to a high convergence rate, irrespective of the nature of boundary conditions [ 222 . 

For smooth functions, pseudospectral methods are much more accurate numerical dif- 
ferentiators than the conventional second- and fourth-order finite-difference approxi- 
mations. However, this superiority of the spectral methods may not be very pro- 
nounced when turbulent flows which often involve small-scale fluctuations are calcu- 
lated. In numerical simulation of two-dimensional, Navier-Stokes equations in a 
periodic box, Herring et al. [232 have systematically compared their results obtained 
with second-order finite-difference and spectral methods. They showed that spectral 
calculations are approximately equivalent in accuracy to finite-difference calcula- 
tions with only twice the resolution in each space dimension. It is quite likely 
that a more favorable comparison can be obtained if fourth-order finite-difference 
methods are used. Another important conclusion from their study was that the accu- 
racy or inaccuracy of spectral methods can be deduced from the computed energy 
spectra, whereas, the spectra obtained from the corresponding finite-difference calcu- 
lations tend to hide their inaccuracy. To illustrate this property of spectral 
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methods, we shall consider two numerical simulations of turbulent channel flow, one 
inadequately resolved and the other with sufficient grid resolution. 

Figure 2a shows the one-dimensional lateral energy spectra in the vicinity of 
the wall = 0.025) from a turbulent channel-flow simulation at Re = 13800, 

Ref. [3]. In this calculation, the pseudospectral method with 128 grid points was 
used in the lateral direction. However, for this Reynolds number the resulting com- 
putational resolution was not adequate to resolve the wall-layer streaks at their 
proper scale. The energy accumulation at the high wave-number portions of 
E^(k 3 ,y ^/6 = 0.025) signals this inadequacy. Note that the longitudinal spectra 
obtained at the same vertical location (Fig. 2b) do not have energy buildup at high 
wave numbers. In this calculation, 64 grid points were used in the longitudinal 
direction. These appear to be sufficient to resolve the streamwise variation of 
turbulent structures (which, incidentally, suggests that the streamwise grid resolu- 
tion estimate given in section 3 may be too stringent (see below)). Figure 2c shows 
the one-dimensional, lateral energy spectra E^(k 3 ,y ^/6 = 0.389) at a distance away 
from the wall. In this region, the finely spaced near-wall structures are absent, 
and no resolution problems are expected. This is reflected in the behavior of 
and the absence of excessive energy buildup at high values of k 3 . Figure 3 shows 

the corresponding near-wall (y /6 = 0.025) one-dimensional spectra from a channel 

w 

flow simulation at the relatively low Reynolds number Re = 3850. In this calcula- 
tion, 128 grid points were also used in the spanwise direction and were apparently 
sufficient to resolve the wall-layer structures. For this case, no excessive energy 
accumulation is evident at the high-wave-number end of E^(k 3 , y^/ ^ = 0.025). One 
should be cautious in using energy spectra as the sole indicator of grid resolution 
adequacy or inadequacy. Insufficient computational resolution may totally suppress 
certain instability mechanisms and the subsequent formation and growth of the corres- 
ponding turbulent structures. This phenomenon may be concealed from energy spectra. 

If the spatial grid resolution is sufficient to resolve all the important scales 
of motion, spectral or pseudospectral methods certainly are the best possible numeri- 
cal differentiators. However, in some cases the accuracy afforded by spectral 
methods should be balanced against inherent inefficiencies in the data-management 
process and difficulties encountered with application of these methods to complex 
geometries. 

4.2 Explicit Time Advancement 

The three momentum Eqs. (4) must be integrated in time, subject to the incom- 
pressibility constraint. Explicit methods offer the advantages of low cost per step 
and ease of formulation and computer programming. In calculations where the wall- 
layer dynamics have been excluded D-, 2 ] , only explicit time advancement has been 


8 



Figure 2. One-dimensional energy spectra at Re 13800. — Ej spectrum 

of the streamwise fluctuating velocity; — — , E 2 spectrum of 

spanwise fluctuating velocity; , E 3 spectrum of normal 

fluctuating velocity, (a) lateral spectra, y^/6 - 0.025; 

(b) longitudinal spectra, - 0.025; (c) lateral spectra, 

y ^^/6 = 0.389. Note that in (a) the ratio of the peak value 
of Ej to its value at the highest resolved wave number is 
only 4.2. 


n 







Figure 3. One-dimensional lateral energy spectra at 

Re = 3850. y^/6 = 0.025 (see caption of Fig. 2). 

used. In these simulations, due to the use of relatively few uniformly or nearly 
uniformly spaced grid points, the stability restriction (especially those of the 
diffusion type) on time step is not severe. Both leapfrog and second-order Adaras- 
Bashforth methods have been employed. The latter method has better overall accuracy 
and stability properties and is more popular now. 

To enforce the incompressibility condition at the next time step, usually, 
the Poisson equation (Ic) for pressure is used rather than the continuity equation. 
For several flow geometries of interest, noniterative, elliptic solvers are available 
C24D for exact solution (to within round-off errors) of the discretized Poisson equa- 
tion. It is important, however, to note that the f inite-dif ferenced operator 

for the Poisson equation cannot be chosen arbitrarily [25 j . In order to ensure com- 
pliance with the incompressibility condition, the numerical gradient operator used 
to approximate 3P/8x^ in Eq. (4) and the divergence operator must be the same. 

(When a staggered grid is used, a combination of forward and backward difference 
schemes for the divergence and the gradient operators also leads to the incompressi- 
bility condition.). Except when second-order finite-difference methods on a stag- 
gered grid are used, there is some ambiguity with the choice of boundary condition 
for the Poisson equation (Ref. [19]). Usually, one uses the Neumann boundary condi- 
tion for pressure (which is obtained from the normal momentum equation) . It is also 
possible to obtain a Dirichlet condition from the boundary evaluation of the tangen- 
tial momentum equations. In general, the Neumann and Dirichlet problems for pressure 
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may not have the same solution. It can be shown [19^ that, when pseudospectral 
methods are used in the direction normal to the boundary, in order to satisfy the 
boundary conditions and the equation of continuity at the wall, it is imperative 
that both Neumann and Dirichlet conditions for pressure be satisfied. However, only 
one of them can be used to solve the Poisson equation. This inconsistency leads to 
serious numerical difficulties if pseudospectral methods are used. With standard 
finite-difference techniques, although the above ambiguity is still present, the 
corresponding numerical difficulties can be avoided [^26, 273. If second-order 
finite-difference methods are used in conjunction with staggered grids, the incom- 
pressibility condition at the interior cell adjacent to the boundary provides the 
additional boundary relation needed to solve the system of algebraic equations for 
pressure. In this case it is not necessary to extract pressure boundary conditions 
from the momentum equations. 

4. 3 Partially Implicit Time Advancement 

The calculations that have extended to the wall and applied the no-slip bound- 
ary conditions have used semi-implicit numerical methods. These numerical schemes 
circumvent the prohibitive time-step restrictions arising from the viscous term and 
the necessity of using very fine mesh spacing in the vicinity of the wall. Moreover, 
when spectral methods are used to approximate the derivatives in the direction normal 
to the boundaries, implicit methods provide the means for convenient imposition of 
boundary conditions. This is because, in contrast to the fully explicit methods, 
at each time step the problem is treated as a boundary value rather than an initial 
value problem. 

The flows considered to date have been restricted to those that are homogeneous 
in two space dimensions. The direction of inhoraogeneity is normal to the wall(s) . 
Specifically, plane channel flow, pipe flow, axial flow between two concentric cyl- 
inders, plane and circular Couette flows have been simulated numerically. For these 
cases, the use of periodic boundary conditions in the homogeneous directions allows 
the application of Fourier transforms, which alleviate the need for split or 
factored- type algorithms. Perhaps, the most direct approach for the solution of 
Eq. (4) and the continuity equation is to solve them simultaneously. For simula- 
tion of turbulent channel flow, Moin and Kim [3, 19j used the second-order Adams- 

Bashforth method for H. and the Crank-Nicolson method for 3P/3x, and 3^u./3x.3x. 

1 113 1 

in Eq. (4). This, together with the continuity equation at the new time step, 
n + 1, led to a system of four coupled, linear, partial-differential equations for 
u^ and P of the form 


- f n+1 T,n+K -f n n-1 _n>. 

L(u^ , P ) = F(u^, u^ , P ) 


(5) 
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where the superscripts denote the time step. Fourier transforming Eq. (5) in the 
xj and X3 directions produces linear, ordinary differential equations of the form 

_2,(-n+l ^ ^n+l^ ^ - 

where denotes Fourier transform. To solve these ordinary differential equations, 
both finite-difference operators on a modified staggered grid C3] and Chebyshev 
polynomial expansions D-91 were used to approximate 3/9x2 3^/3x2. The result 

is a system of algebraic equations for the Fourier transform of the dependent vari- 
ables at the new time step. For the case where finite-difference operators were 
used, this system is of block— tridiagonal form; and for the case where Chebyshev 
polynomials were used, it is nearly block— tridiagonal. Both systems can be solved 
with 0(N2) operations, where N2 is the number of mesh points in the X2-direction. 

For their study of transition to turbulence in plane channels and Couette flows, 
Orszag and Kells [17] used a three-step fractional step method. Chebyshev poly- 
nomials and Fourier decomposition were used to represent the dependent variables 
spatially. In the first step, the Adams-Bashforth method is used for the nonlinear 
terms, H.. The result is a set of intermediate velocity field, . Next, the 

pressure correction (incompressibility condition) is applied, leading to another set 
of intermediate velocity field u^"*"^ that satisfies the continuity equation 

3u^'*’V9x = 0. This step involves solving a Poisson equation for with the 

i i _ 

boundary conditions u" = 0 at the walls. The velocity field at the new time 
step is obtained by applying the viscous correction which involves the solution of 
three Helmholtz equations for the velocity field. The velocity boundary conditions 
are applied at this stage. This method has a global error of order 0 [At^+ (1/Re) At]. 
Thus, strictly speaking, it is a first-order method. In addition, the velocity 
field at the new time step does not satisfy the continuity equation. Only the inter- 
mediate velocity field, fi. , is solenoidal. If the last two steps were interchanged, 
the velocity field would be divergence-free, but the boundary conditions could not 
be enforced. Apparently, having the exact boundary conditions was preferred over 
the continuity equation. Another characteristics of note is that, to solve the 
Helmholtz equations in the final step using Chebyshev-Fourier expansions, one must 
transfer the forcing terms in these equations to the wave (Chebyshev-Fourier) space. 
These terms include . Thus, to carry out the transformation, u^^^ must be 

defined at the boundaries as well as in the interior of domain. It is necessary, 

_ T . • T . . ::n+l 

therefore, to concoct boundary conditions for the intermediate velocities, 

Recently, Leonard [5] developed a partially implicit method based on a special 
vector-function decomposition. An important feature of these vector functions is 
that each vector is solenoidal and satisfies the boundary conditions. IThen this 
series expansion for the velocity vector is used in Eq. (4) and the inner product of 
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the result with a set of adjoint vectors is formed, the pressure terms are elimi- 
nated. In fact, since the velocity bases functions satisfy the continuity equation, 
only two dependent variables per mesh point remain. This leads to considerable 
savings in computer-memory requirements. Like the previous methods, the nonlinear 
terms are treated by an explicit method, whereas, the Crank-Nicolson scheme is used 
for the viscous terms. Leonard and Wray [2S2 have applied this method to the pipe- 
flow problem. They rigorously treat the behavior of the flow variables near the 
computational singularity at the pipe centerline which leads to the use of the 
Jacobi polynomial expansions in the radial direction. Moser et al. C27l applied 
this method to plane and curved channel flow problems. In the direction normal to 
the walls, a particular combination of Chebyshev polynomials was used as bases 
functions. With this choice, for each Fourier mode, the resulting system of alge- 
braic equations was solved with 0(N2) operations. 

4 . 4 Data Management 

Limitations of the high-speed central memory of presently available computers 
and the large number of grid points required for numerical simulation of turbulent 
flows necessitate the use of secondary memory. Generally, the entire data base 
must reside on secondary memory (SM) and only portions of it are successively trans- 
ferred to the core memory (CM) for processing. In order to minimize the data trans- 
fer (I/O) time, an efficient data-management algorithm should be an essential part 
of each computer program. 

The particular choice of data-management technique depends on the numerical 
method and the computer used for the calculations. However, all the algorithms have 
the objective of overlapping the data transfer from SM to CM with arithmetic opera- 
tions. It is also important to minimize the number of passes over the data base. 

In general, one or two passes are required at each time step. 

As an example of a data-management algorithm, we consider the scheme employed 
in Ref. [3]. In order to solve Eq. (6), a two-pass, double-buffer, data-management 
algorithm was used. In the first pass over the data base, the required pressure- 
velocity data from previous time steps are transferred to CM, and F (Eq. (6)) is 
computed and tTansf erred to SM. Since second— order finite— difference formulas were 
used to approximate the derivatives in the x2~direction, each time only three 
(xi - X3) planes of data are required to compute F in one plane. In the second 
pass, (x2 “ X3) planes of F were transferred to CM, and the block-tridiagonal 
system was solved. Note that, with this algorithm, the data were accessed in two 
different ways — first in the xj - X3 planes and then in the X2 - X3 planes. 
Since the data are stored sequentially in, say, the x^ - X3 planes, the second 
access (x2 - X3 planes) is nonsequential; therefore, it is not as efficient. 
Alternatively, another data-management scheme can be used for solving the same set 
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of equations. With this algorithm, in the first pass, in addition to computing F, 
the forward sweep portion of the Gauss-elimination process can also be performed. 

In the second pass, the backward sweep portion is performed again, using - X3 
planes. Thus, the data base is accessed sequentially. In fact, the entire process 
can usually be accomplished with only one pass through the data base per time step 
C30j . It should be pointed out that, if pseudospectral rather than finite- 
difference methods are used in the X2-direction , two passes through the data base 
(one of them nonsequential) are required. 

We emphasize that the numerical methods discussed here have been largely 
designed for and applied to the simplest wall-bounded, turbulent shear flows, 
namely, those that are homogeneous in two spatial directions. The flow homogen- 
eity has led to the use of periodic boundary conditions as a reasonable approxima- 
tion to the unknown flow conditions at the "open** boundaries. The use of 
Fourier transforms converts the task of solving a partial differential equation to 
that of solving a set of uncoupled ordinary differential equations. This is a very 
significant fringe benefit associated with periodic boundary conditions. In order 
to calculate flows with two or more directions of inhomogeneity, in addition to 
having to specify the unknown (turbulent) inflow and outflow conditions, one should 
use split- or factored-type algorithms. An adaptation of the fractional step method 
appears to be an attractive candidate for this purpose, 

5. Results 

To illustrate the versatility and usefulness of the aforementioned calculations, 
we briefly present here a selection of recent results from the large-eddy simulation 
of wall-bounded turbulent flows. The main calculations were performed on the ILLIAC 
IV computer with 64x63x128 grid points in the x, y, and z directions, respec- 
tively. The total computational time ranged from 20 hours to 92 hours. 

Figure 4 shows the mean velocity profile, (u), from numerical simulation of 
turbulent channel flow at Re = 13800, Ref. 3. Different symbols represent calcula- 
tions with different grid resolutions and different sizes of the computational box 
in the x- and z-directions, where periodic boundary conditions are used. The cal- 
culations have predicted the logarithmic region with the proper slope. The agree- 
ment with the experimental data [3lD is good. The distributions of the Reynolds 
stresses (Fig. 5) and higher-order statistics are also in good agreement with mea- 
surements. The contribution of subgrid scale turbulence to second- and higher-order 
statistical correlations is appreciable only in the vicinity of the walls. This is 
a consequence of the grid-resolution inadequacy in this region to represent the wall- 
layer turbulence structures at their proper scale. Figure 6 shows a contour plot of 
the instantaneous normal component of vorticity fluctuations, ujo = (3w/3x - 3u/9z) 
in an x - z plane close to the wall (y"*" ^ 6). It is clear that, in accordance 
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Figure 4. Mean velocity profile from four computed cases 
and comparison with experimental data. Four 
results are obtained from using different com- 
putational grids. 



O COMPUTATION 
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COMTE-BELLOT 

(1963), Rt- 64000 
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Fig. 5 Resolvable turbulence intensities and comparison 
with experimental data. The intensities are non- 
dimensionalized with the shear velocity, u^ . 
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Figure 6* Contours of cj 2 in the (x, z)-plane at = 6.26. The rectangle on the 
lower right-hand corner of the figure represents the computational grid 
cell in the (x, z) planes. The streamwise extent of the figure is 2tt6 
and its spanwise extent is 7 t6, 

with experimental measurements [6l , this region is composed of flow structures which 
are long in the flow direction and narrow in the spanwise direction. The rapid 
spanwise variation of cj 2 is due to the existence of elongated regions of high- 
speed fluid ((u - (u)) > 0) located adjacent to the low-speed regions [6, 3]. 

This figure is a vivid display of that particular characteristic of wall-bounded 
turbulent shear flows which requires a large number of grid points in the lateral 
direction. As was pointed out in section 3, in this calculation, the spanwise grid 
resolution was not adequate to resolve the wall-layer streaks at their proper scale. 
However, it is quite significant that, in spite of this, the computed flow field 
did display the streaky structures but at a larger scale. 

The data generated from these calculations are currently being used to study 
the physical structure and dynamics of turbulent channel flow. In one study, for 
example, our aim is to identify large-scale, energetic structures in the flow field. 
In particular, we wanted to investigate the frequency and dominance of horseshoe or 
hairpin vortices that have been observed to originate at the wall and extend to 
outer regions with the characteristic inclination angle 40° - 50° (see, e.g.. Ref. 
C32j.) The vorticity field at several points in time was computed. At each grid 
point in various x - z planes, the angle 0 = tan ^(cl) 2 /o)i), and 
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I“l2l “ ^ 2 ^ was calculated. Figure 7 shows the resultant distribution at 

= 0.2. The contribution from each grid point was weighted by |t^l2l- Indeed, 
the distribution attains its maximum at 0 = 45®. However, the probability of find- 
ing vorticity vectors with inclination angle in the range 0 < 0 < 90° is also 
appreciable. A detailed description of the results of this study will be presented 
elsewhere. In another investigation, Kim CSSD has used the computed velocity- 
pressure field to examine the structure of the flow by conditional sampling tech- 
niques. Figure 8 shows the conditionally averaged, streamwise velocity obtained by 
using the VITA technique C7j . It is remarkably similar to the experimental results 
of Blackwelder and Kaplan C73 • The figure displays the burst and sweep events and 
their vertical extent. Kim has extended the experimental findings utilizing condi- 
tionally averaged pressure, vorticity, and the spanwise velocity component. 

By simple modifications of the channel flow code described above, the effects 
of transpiration and spanwise rotation on the flow were computed. In the former case, 
uniform blowing through one wall and uniform suction at the same rate was applied 
through the other wall. Figure 9 shows that, in agreement with experimental measure- 
ments, the calculations predict the wall-shear-stress diminution caused by blowing 



d, deg 

Figure 7. Distribution of the inclination 

angle of vorticity vectors (weighted 

by |wi 2 |) at y /5 = 0.2. 
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Figure 8. Conditionally averaged streamwise 

velocity fluctuations at various dis- 
tances from the wall. 

and its augmentation resulting from suction. Other features of the computed flow 
field, such as the effect of transpiration on the distribution of Reynolds stresses, 
are also in agreement with measurements. In the case of the rotating channel, the 
computational results CSSD are in good agreement with the experimental data C36] and 
reproduce the detailed structural features of the flow as observed by flow- 
visualization techniques. The primary effect of rotation is that the flow is stabi- 
lized on one wall and destabilized on the other. The skin friction is reduced on 
the ’’stable wall” and is increased on the ’’unstable wall” (Fig. 10). One of the 
objectives of these studies was to examine the relationship between the changes in 
skin friction and flow structure in the vicinity of the wall. In both flows, it 
has been shown that there is a definite correlation between the characteristic 
dimensions of the wall-layer streaks and viscous drag. 

6. Concluding Remarks 

The results of the numerical simulation of wall-bounded turbulent shear flows 
have been most encouraging. For geometrically simple cases, it has been possible to 
predict many of the statistical and time-dependent features of the flows considered. 
The potential of these calculations for increasing our understanding of the physics 
of turbulent boundary layers is just beginning to be tapped. 
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Figure 9* Ratio of friction coefficient to friction co- 
efficient without transpiration, plotted as a 
function of a modified blowing parameter hf [34]. 

O , computation; correlation curve of the 

experimental data; , linear extrapolation of 

the experimental correlation curve. 

The large number of grid points required to resolve the wall-layer structures 
is a formidable hindrance to numerical simulation of high Reynolds number flows. 

The grid-imbedding technique can significantly ease this burden. Consideration of 
the bursting frequency and the frequency spectra of streamwise velocity fluctua- 
tions [37, 13] indicates that, for moderate Reynolds numbers, it may be possible to 
perform accurate LES calculations with larger time steps than is currently permitted 
by numerical stability restrictions. Thus, further improvements in numerical 
methods are likely to yield high dividends. 

Extension to more complex geometries requires progress in our ability to pre- 
scribe turbulent inflow and outflow boundary conditions. Another practical diffi- 
culty in calculating these inhomogeneous flows is acquiring adequate ensemble 
averages. In contrast to homogeneous flows, one cannot integrate turbulence quanti- 
ties over spatial grid points to secure a better statistical sample. To obtain 
adequate statistics, the only resort is time-averaging or, in the case of 
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Figure 10. Ratio of the shear velocity to the shear 
velocity without rotation, plotted as a 
function of the rotation number. A , com- 
putation p5j ; , curve fit to the 

experimental data [36j . 


nonstationary flows, averaging over several independent calculations. In both 
cases, the computational cost is significantly higher than for flows with one or two 
directions of homogeneity. 
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